//
//  homotopy.cpp
//  EMT_V7_LSAP_Guess
//
//  Created by Duong Ngo on 7/4/17.
//  Copyright © 2017 Duong Ngo. All rights reserved.
//

#include "homotopy.hpp"
#include "foc_problem.hpp"
#include <fstream>
#include <iomanip>


void store_output_homotopy (const std::vector<double>& result, const double& chi, const std::string& file_name)
{
    std::ofstream file;
    file.open(file_name.c_str());
    
    file << T*var_foc << " " <<  chi << std::endl;
    for (int t=0; t<=T; ++t){
        for (int i=0; i<var_foc; ++i){
            file << t << " " << i << " " << std::setprecision(25)<< result[t*var_foc+i] << std::endl;
        }
    }
    file.close();
}

//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------

void conduct_homotopy (const std::vector<double>& kappa_start, const std::vector<double>& kappa_end, const std::vector<double>& x_start, const std::vector<double>& x_end,  const std::vector<double>& rn_start, const std::vector<double>& rn_end, const std::vector<double>& ss, const std::string& file_name, std::vector<double>& x_init)
{
    int N_iter = 100;
    int dem =0;
    double chi=0; //Homotopy coefficient
    double acceptable_error = 1e-6;
    
    while (dem <= N_iter && max_error < acceptable_error){
        printf("*********************The iteration: %2u *********************\n", dem);
        chi= double(dem)/N_iter;
        for (int t=0; t<T+1; ++t){
            vec_kappa[t]= chi*kappa_end[t]+ (1-chi)*kappa_start[t];
            vec_x[t]= chi*x_end[t]+ (1-chi)*x_start[t];
            vec_rn[t]= chi*rn_end[t]+ (1-chi)*rn_start[t];
        }
        foc_solve(x_init, ss);
        if (max_error < acceptable_error){
            x_init = result;
            store_output_homotopy(result, chi, file_name);
        }
        else {
            for (int j=1; j<= 100; ++j){
                chi = double(dem-1)/N_iter+ double(j)/100/N_iter;
                for (int t=0; t<T+1; ++t){
                    vec_kappa[t]= chi*kappa_end[t]+ (1-chi)*kappa_start[t];
                    vec_x[t]= chi*x_end[t]+ (1-chi)*x_start[t];
                    vec_rn[t]= chi*rn_end[t]+ (1-chi)*rn_start[t];
                }
                foc_solve(x_init, ss);
                if (max_error < acceptable_error){
                    x_init = result;
                    store_output_homotopy(result, chi, file_name);
                }
                else {
                    break;
                }
            }
        }
        ++dem;
    }

}

//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------

/** Homotopy from steady state to homo_kappa*/
void homotopy_from_ss_to_homo_kappa (const std::vector<double>& ss, std::vector<double>& kappa_start, std::vector<double>& kappa_end, std::vector<double>& x_start, std::vector<double>& x_end,  std::vector<double>& rn_start, std::vector<double>& rn_end,  std::vector<double>& x_init)
{
    for (int t=0; t<T+1; ++t){
        kappa_start[t]= kappa_homotopy;
    }
    kappa_end[0] = kappa_shock;
    for (int t=1; t<T+1; ++t){
        kappa_end[t]= rho*kappa_end[t-1]+ (1-rho)*kappa_homotopy;
    }
    
    for (int t=0; t<T+1; ++t){
        x_start[t] = 0;
    }
    x_end[0]= 0; x_end[1]= x_shock;
    for (int t=2; t<T+1; ++t){
        x_end[t]= rho_x*x_end[t-1];
    }
    
    for (int t=0; t<T+1; ++t){
        rn_start[t]= Rn;
        rn_end[t]= Rn;
    }
    
    std::string file_name ("//Users//duongngo//Documents//Research//E-Monetary Theory//Paper_and_code//EMT_V6//C code//Output//EMT_ms_rule_homo_kappa.txt");
    
    conduct_homotopy(kappa_start, kappa_end, x_start, x_end, rn_start, rn_end, ss, file_name, x_init);

}

//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------

/** Homotopy from homo_kappa to kappa */
void homotopy_from_homo_kappa_to_kappa (const std::vector<double>& ss, std::vector<double>& kappa_start, std::vector<double>& kappa_end, std::vector<double>& x_start, std::vector<double>& x_end,  std::vector<double>& rn_start, std::vector<double>& rn_end,  std::vector<double>& x_init)
{
    kappa_start[0]= kappa_shock;
    kappa_end[0] = kappa_shock;
    for (int t=1; t<T+1; ++t){
        kappa_start[t]= rho*kappa_start[t-1]+ (1-rho)*kappa_homotopy;
        kappa_end[t]= rho*kappa_end[t-1]+ (1-rho)*kappa;
    }
    
    x_end[0]= 0; x_end[1]= x_shock;
    for (int t=2; t<T+1; ++t){
        x_end[t]= rho_x*x_end[t-1];
    }
    x_start= x_end;
    
    for (int t=0; t<T+1; ++t){
        rn_start[t]= Rn;
        rn_end[t]= Rn;
    }
    
    std::string file_name ("//Users//duongngo//Documents//Research//E-Monetary Theory//Paper_and_code//EMT_V6//C code//Output//EMT_ms_rule_without_rn.txt");
    
    conduct_homotopy(kappa_start, kappa_end, x_start, x_end, rn_start, rn_end, ss, file_name, x_init);
}

//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------

/** After "T_zero" periods, the central bank raise "Rn" to 1./betab */
void homotopy_for_vec_rn (const int& T_zero, const std::vector<double>& ss, std::vector<double>& kappa_start, std::vector<double>& kappa_end, std::vector<double>& x_start, std::vector<double>& x_end,  std::vector<double>& rn_start, std::vector<double>& rn_end,  std::vector<double>& x_init)
{
    kappa_end[0]= kappa_shock;
    for (int t=1; t<T+1; ++t){
        kappa_end[t]= rho*kappa_end[t-1]+ (1-rho)*kappa;
    }
    kappa_start = kappa_end;
    
    x_end[0]= 0; x_end[1]= x_shock;
    for (int t=2; t<T+1; ++t){
        x_end[t]= rho_x*x_end[t-1];
    }
    x_start= x_end;
    
    for (int t=0; t<T+1; ++t){
        rn_start[t]= Rn;
        rn_end[t]= Rn;
        if (t> T_zero && t< T_rn_back){
            rn_end[t]= 1./beta_b;
        }
    }
    
    std::string file_name ("//Users//duongngo//Documents//Research//E-Monetary Theory//Paper_and_code//EMT_V6//C code//Output//EMT_ms_rule_raise_rn.txt");
    
    conduct_homotopy(kappa_start, kappa_end, x_start, x_end, rn_start, rn_end, ss, file_name, x_init);
}






